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Abstract 

A simple 4-node assumed-stress hybrid quadrilateral shell element with rotational or 
“drilling” degrees of freedom is formulated. The element formulation is based directly 
on a 4-node element. This direct formulation requires fewer computations than a simi- 
lar element that is derived from an “internal” 8-node isoparametric element in which the 
midside degrees of freedom are eliminated in favor of rotational degrees of freedom at the 
corner nodes. The formulation is based on the principle of minimum complementary en- 
ergy. The membrane part of the element has 12 degrees of freedom including rotationa 
degrees of freedom. The bending part of the element also has 12 degrees of freedom. The 
bending part of the element uses the Reissner-Mindlin plate theory which takes into ac- 
count the transverse shear effects. Quadratic variations for both in-plane and out-of-plane 
displacement fields and linear variations for both in-plane and out-of-plane rotation fields 
are assumed along the edges of the element. The element Cartesian-coordinate system is 
chosen such as to make the stress field invariant with respect to node numbering. 1 e 
membrane part of the stress field is based on a 9-parameter equilibrating stress field, while 
the bending part is based on a 13-parameter equilibrating stress field. The element passes 
the patch test, is nearly insensitive to mesh distortion, does not “lock,” possesses the de- 
sirable invariance properties, has no spurious modes, and produces accurate and reliable 

results. 


Introduction 

Finite element researchers face what seems to be an endless challenge to formulate simple 
3-node and 4-node shell elements that are free from the usual deficiencies, such as locking, 
sensitivity to mesh distortion, non-invariance, and spurious modes. From the inception of 
the standard 4-node isoparametric element, researchers realized that this element exhib- 
ited severe locking and was very sensitive to mesh distortion. Ever since, researchers have 
considered a variety of methods to overcome these deficiencies. These methods have elim- 
inated some of the shortcomings of the standard 4-node isoparametric element. However, 
some new difficulties such as non-invariance and spurious modes were introduced, borne 
of the milestones in the quest for a defect-free 4-node element are: 

(1) Assumed-stress hybrid elements (Pian[l]). 

(2) Reduced integration (Zienkiewicz et al.[2] and Pawsey and Clough[3]). 

(3) Incompatible elements (Wilson et al.[4] and Taylor et al.[5]). 


Another method of attacking the shortcomings of membrane elements is to include the 
nodal rotational or “drilling” degrees of freedom in the element formulation. In early 
attempts, these rotational degrees of freedom were used in cubic displacement functions. 
However, Irons and Ahmad demonstrated that this approach had serious deficiencies [6]. 
The elements formed in this manner force the shearing strain to be zero at the nodes, and 
these elements do not pass the patch test, which could produce erroneous results in some 
structural analysis problems. Recently researchers have used these rotational degrees of 
freedom in quadratic displacement functions with more success[7-ll]. In previous papers, 
this latter method has been employed in the following way. First, the element is inter- 
nally assumed to be an 8-node isoparametric element with 4 corner nodes and 4 midside 
nodes each having two displacement degrees of freedom, and the stiffness matrix associ- 
ated with this “internal” element is calculated. Then, this stiffness matrix is condensed to 
that corresponding to a 4-node element with 12 degrees of freedom by associating the dis- 
placement degrees of freedom at the midside nodes with the displacement and rotational 
degrees of freedom at the corner nodes. MacNeal[9] has used this approach to develop 
a 4-node displacement-based membrane element with selective reduced-order integration. 
Yunus et al.[10] have also used this method to develop an assumed-stress hybrid/mixed 
membrane element. Aminpour[ll] has also used this method to develop an assumed-stress 
hybrid/mixed shell element. 

In this paper, a 4-node assumed-stress hybrid quadrilateral shell element with rotational 
degrees of freedom is presented. The formulation is based directly on a 4-node element 
from the beginning in contrast to elements whose formulations began with an “internal” 
8-node element. Formulating the element in this manner bypasses the formation of the 
stiffness matrix for an 8-node isoparametric element and the subsequent transformation of 
this stiffness matrix to that corresponding to stiffness matrix of a 4-node element. This 
method is advantageous in that the element formulation is more direct and savings in 
computations are accrued. Results are presented for several standard test problems to 
establish the robustness of this element. 


Hybrid Variational Principle 

The classical assumed-stress hybrid formulation of Pian[l] is based on the principle of 
minimum complementary energy. The displacements are described on the element bound- 
ary and an equilibrating stress field is described over the the domain of the element. It 
was later recognized that the same method may be derived from the Hellinger-Reissner 
pnnciple[12-14]. However, in the Hellinger-Reissner principle, the stress field does not have 
to satisfy the equilibrium equations a priori, and the displacement field has to be described 
over the domain of the element and not just on the boundaries. The stress field would 
then satisfy the equilibrium equations only in a variational sense. Therefore, the stress 
field may be described in the natural-coordinate system of the element which would make 
the element less sensitive to mesh distortion, and a proper selection of the stress field would 
make the element invariant with respect to node numbering. Because of these desirable 
properties, researchers have developed assumed-stress hybrid/mixed elements using the 
Hellinger-Reissner principle. For example, the membrane element in reference [10] and the 
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shell element in reference [11] were both developed using the Hellinger-Reissner principle. 
However, an assumed-stress hybrid 4-node shell element similar to that of reference [11] 
may also be easily formulated using the minimum complementary energy principle with 
the advantage being that only the displacements on the boundary of the element enter 
into the formulation. As such, the formulation is then based directly on a 4-node element 
rather than internally formulated as an 8-node element and then condensed to a 4-node 
element. 

The invariant properties of the element is preserved by proper choice of a local element 
Cartesian-coordinate system. The local element Cartesian-coordinate system is shown in 
Figure 1 and is obtained by bisecting the angles formed by the diagonals of the element. 
The axes of this coordinate system are approximately parallel to the edges of the element 
for non-rectangular geometries (e.g., tapered and skewed elements) which would make the 
element less sensitive to mesh distortion. Upon node renumbering, this coordinate system 
is rotated by 90° increments. For example, if the connectivity for the element shown in 
Figure 1 is changed from 1-2-3-4 to 2-3-4-1, then the element local x-y axes in Figure 1 
would rotate by 90°. Therefore, selecting stress fields that are invariant with respect to a 
90° rotation would make the element invariant with respect to node numbering. 

The formulation of the element presented herein is based on the principle of minimum 
complementary energy. The details of the assumed-stress hybrid formulation using the 
minimum complementary energy principle have been extensively discussed in the literature, 
(e.g., see reference [1]), and hence, only a brief outline is given herein for completeness. 
The variational functional is given as 

n = — ^ f cr T DcrdV + f a T nudS — f u T t 0 dS (1) 

2 Jv Js Js t 

where D is the compliance matrix of the material, a is the stress array, u is the displacement 
array, to is the prescribed traction array, matrix n consists of the components of the 
outward unit normal to the boundary of the element such that n T er = t (traction array), 
V is the domain of the element, S is the boundary of the element, and S t is the part of 
S where to is specified. The assumed-stress hybrid formulation is based on assuming an 
equilibrating stress field in the interior of the element as 

* = P/3 (2) 

and assuming the displacement field only on the boundary of the element as 

u — Nq (3) 

where the matrices P and N consist of the appropriate interpolating functions for stresses 
and displacements, respectively, and the coefficients /3 and q are the unknown stress pa- 
rameters and nodal displacements and rotations, respectively. 

The expressions for stresses in equation (2) and displacements in equation (3) are substi- 
tuted into the functional n of equation (1) and the variation of the functional with respect 
to the internal unknowns (3 is set to zero. This stationary condition gives 

/3 = H -1 Tq (4) 
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where 


H = / P T DPdF 
Jv 

(5) 

T = J P T nN dS 

(6) 

Upon substitution of the expression for (3 in equation (4) into the functional II of equation 
(1) and a subsequent variation of the functional with respect to the nodal unknowns q 
yields 

Kq = F (7) 

where the stiffness matrix K is given by 

K = T t H -1 T 

(8) 

and the generalized force array F by 


F = / N T t 0 d5 
Js* 

(9) 


Formulation of 4-Node Quadrilateral Element 

Displacement Field Description 

As discussed previously, all three displacement components on the element boundary are 
assumed to vary quadratically and all three rotational components to vary linearly. The 
expressions for these boundary displacements and rotations were derived in detail in refer- 
ence [11] and only the fined results are given herein. The in- plane boundary displacements 
on edge 1 of Figure 1 are given by 

« = 5(1 - {)«> + 5(1 + 0<*> + - ( 2 )(t>x - «„) 

? ? A l (10) 

» = j(l - f)»l + j(l + - -ji(l - - «,.) 

and the out-of-plane boundary displacement and rotations on edge 1 of Figure 1 are given 
by 

1 1 A 7 / At* 

» = j(i - {)«>. + 2<i + - x (1 " £2)(9 * 2 ~ hi) + t (1 - < - *»•) 

— 0^*1 + 2 + 0^*2 ( 11 ) 
= ^(1 — + ^(^ ^)^» 2 

where, Axi and A y\ are the Ax and A y of edge 1 with respect to the reference local 
element x-y coordinate system (e.g., Axj = x 2 — Xi) and £ is a non-dimensional coordinate 
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on edge 1 such that £ = —1 at node 1 and £ = +1 at node 2. It is worth mentioning here 
that the true nodal normal rotations are, of course, given by — §^) evaluated at the 

nodes. Hence, the terms 0 Z ; are not true nodal rotations, and they may be referred to 
as “rotational connectors” [7]. The description for the displacements and rotations on the 
other edges of the element are readily obtained. 

The description for the in-plane rotation 8 Z is similar to the out-of-plane rotations 9 X and 
9y. However, the description for the in-plane rotation 0 Z is not shown here and it does not 
enter into the membrane formulation, while both out-of-plane rotations 8 X and 9 y enter into 
the bending formulation. Therefore, both membrane part and bending part of the element 
are formulated in the same manner. All three displacement components are quadratic 
functions, while all three rotations are linear functions. This conformity in the order of the 
approximating polynomials for all displacement components and all rotational components 
is very desirable in analysis of shell problems. The displacement and rotation descriptions 
in equations (10) and (11) allow for in-plane shearing strain and transverse shearing strains, 
respectively. This feature is in contrast to cubic interpolations of in-plane or out-of-plane 
displacements which force the in-plane shearing strain or the transverse shearing strains 
to be zero at the element nodes. As discussed previously, elements constructed using 
quadratic interpolation have been more successful. These elements pass the patch test 
which is a necessary condition for convergence to the correct solution. The elements using 
cubic interpolation, on the other hand, do not pass the patch test and perform poorly for 
some structural analysis problems [6]. 

Equations (10), when extended to all four sides of the element, indicate that for the 
membrane part of the element 4 in-plane rotational degrees of freedom in addition to 
the usual 8 in-plane displacement degrees of freedom are required to express the in-plane 
displacements as quadratic functions. The bending part of the element, on the other hand, 
is formulated in terms of the usual 12 out-of-plane displacement and rotational degrees 
of freedom and no additional degrees of freedom are required to express the out-of-plane 
displacement as a quadratic function. 

Therefore, the membrane part of the element has 12 degrees of freedom. Two of these are 
the in-plane translational rigid body motions and one is the in-plane rotational rigid body 
motion. Of the nine remaining degrees of freedom, three represent the constant strain 
states, five represent higher-order strain states, and the final degree of freedom represents 
a special type of zero-energy or “spurious” mode. As discussed in reference [11], this zero- 
energy mode is of a special type and is different from other spurious mechanisms such 
as the hour-glass mode. This zero-energy mode is associated with a state of zero nodal 
displacements and equal nodal rotations 0 2 , which renders the in-plane displacements 
u and v in equation (10) to be zero on all edges of the element. This mode is shown 
in Figure 2 using a cubic interpolation for displacements and may be called the “zero 
displacement” mode[7]. As discussed in reference [11], this zero-energy mode appears 
because the displacements are based on the differences of the nodal rotations and not the 
nodal rotations themselves. Therefore, the membrane part of the element has, in fact, 
only 3 independent rotational degrees of freedom but is expressed in terms of 4 rotational 
degrees of freedoms. Hence, one of the rotational degrees of freedom is superfluous and 
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must be eliminated. This zero-energy mode may be eliminated simply by prescribing at 
least one rotational degree of freedom in the entire finite element model of the structure. 

As discussed earlier, the bending part of the element also has 12 degrees of freedom. One 
of these is the out-of-plane translational rigid body motion and two of these are the out- 
of-plane rotational rigid body motions. Of the nine remaining degrees of freedom, three 
represent the constant curvature states, two represent the constant transverse shear strain 
states, and the other four represent higher-order strain states. No zero-energy modes 
are associated with the bending part of the element despite the fact the out-of-plane 
displacement is also expressed in terms of the differences of nodal rotations. The reason is 
that the out-of-plane rotations 9 X and 9 y enter into the bending formulation to account for 
each one of the out-of-plane nodal rotations, but the in-plane rotation 0 Z does not enter 
into the membrane formulation. 

Stress Field Description 

The stress field should be selected in such a manner that no spurious zero-energy mode 
is produced. A spurious zero-energy mode is produced when the product of a selected 
stress term and the strains that are derived from the displacement functions produces 
zero strain energy under a particular, but not trivial, deformational displacement field. In 
order to avoid spurious zero-energy modes, each independent stress term must suppress 
one independent deformation mode. Therefore, the minimum number of stress terms 
required is equal to the number of degrees of freedom of the element less the number 
of rigid body modes. Spurious zero-energy modes generally occur for regular geometries 
such as rectangular planar elements and brick solid elements and disappear for irregular 
geometries [15]. As discussed previously, the stress fields (membrane and bending) are 
expressed in a proper Cartesian-coordinate system and selected in such manner as to 
remain invariant upon node renumbering. This coordinate system is shown in Figure 1 
and the stress fields that are expressed in this coordinate system must remain invariant 
under a 90° rotation to remain invariant under node renumbering. The selected stress 
fields must also satisfy the equilibrium equations in order to be used in the functional II 
of equation (1). 

As discussed previously, the membrane part of the element has 12 degrees of freedom, three 
of which are due to the in-plane rigid body motions. Therefore, a stress field with a min- 
imum of 9 independent parameters is needed to describe the membrane stress (resultant) 
field. The following equilibrating stress (resultant) field is considered for the membrane 
part 

N x = 0i + 042/ + 06* + 08l/ 2 

N y = 02 + 05* + 072/ + 0 9* 2 ( 12 ) 

N X y = 03 - 062/ - 07* 

This stress (resultant) field is expressed in the local element Cartesian-coordinate system 
shown in Figure 1 and is similar to that proposed by reference [11]. However, in refer- 
ence [11] the Hellinger-Reissner principle was used, and the stresses were expressed in the 
natural-coordinate system. The first five terms of the stress field in equation (12) represent 
the stress field that was used in the original 4-node (see reference [1]) assumed-stress hybrid 
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membrane element with 8 degrees of freedom which did not include any normal rotational 
degrees of freedom. The remaining four terms are present to suppress the four rotational 
degrees of freedom present in this formulation. As discussed in reference [11], this se- 
lection of stress field produces no spurious zero-energy modes for the assumed in-plane 
displacement field described in equation (10). 

As discussed earlier, the bending part of the element has 12 degrees of freedom, three 
of which are due to the out-of-plane rigid body motions. Therefore, a stress field with a 
minimum of 9 independent parameters is needed to describe the bending stress field. The 
following equilibrating stress (resultant) field is selected here for the bending part 

M x = P i + Piy + P&x + P&xy 


My = p 2 + Psx + p 7 y + p 9 xy 

M xy = p 3 + Pi ox + p n y + ^/W 2 + \pisy 2 


(13) 


The transverse shearing forces Q x and Q y are obtained using the equilibrium equations 


Q* = 


dM x dM. 


dx 


+ 


xy 


_ PMxy 

Wy ~ dx 


+ 


dy 

dMy 


(14) 


dy 


(15) 


which gives _ _ 

Qx = (Pe + Pn) + (Pa + Pi*)y 
Qy = {Pf + Pio) + (Pa + Pi2)x 

This stress (resultant) field is expressed in the local element Cartesian-coordinate system 
shown in Figure 1 and is similar to that of reference [11]. However, in reference [11], 
the Hellinger-Reissner principle was used and the stresses were expressed in the natural- 
coordinate system. The stress (resultant) field given by equations (13) and (15) is obtained 
by integrating, through the thickness, the stress field that was derived in references [16-18]. 
The stress field in references [16-18] was derived by expressing the stress components as 
power series in the plate thickness, substituting these stresses into the continuum equations 
of elasticity, and equating the coefficients of like powers of the plate thickness. This 13 
parameter selection of stresses for the bending part is less sensitive to geometric distortion 
than a 9 parameter selection obtained from a degenerate solid model[19]. This selection of 
stresses produces no spurious zero-energy modes. It is observed that both the membrane 
and bending stress (resultant) fields remain invariant upon node renumbering. 

In this paper, the Reissner-Mindlin plate theory is used for the bending part. The bending 
part of the element is of class C° and takes into account the effects of transverse shear 
deformations by assuming constant transverse shear strains through the thickness of the 
plate. This assumption means that the transverse shear stresses are also constant through 
the thickness of the plate. However, generally the transverse shear stresses are zero on 
the plate surfaces. Therefore, a parabolic variation of transverse shear stresses and strains 
through the plate thickness is more reasonable. To account for this discrepancy, a static 
correction factor of 5/6 is included in the transverse shear strain energy, see Reissner[20]. 
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Other Element Matrices 

Formulating the element in the manner described above, provides a direct derivation of 
the linear stiffness matrix for a 4-node assumed-stress hybrid shell element. However, this 
methodology is not carried over to other element matrices because, for element matrices 
such as mass and geometric stiffness matrices the displacement functions over the domain 
of the element are required. The derivation of these element matrices are based on the 
“internal” 8-node element described in reference [11]. 

The surface traction and pressure loads are assumed to vary bilinearly over the element 
surface and the force vector F is calculated for the “internal” 8-node element. This force 
vector is then condensed to that of the 4-node element using the approximations for the 
midside degrees of freedom described in reference [11]. The calculation of the force vector 
F for line loads, however, are conducted based directly on the 4-node element by assuming 
linear variation of line loads along the edges of the element and using the displacement 
and rotation shape functions of equations (10) and (11). 

Numerical Results 

The performance of the 4-node quadrilateral shell element developed in this paper is eval- 
uated in this section. The element has been implemented in the NASA Langley CSM 
Testbed software system[21] using the generic element processor template[22]. Selected 
test problems are reported including the patch test, the straight cantilever beam, the 
curved cantilever beam, Cook’s tapered and swept panel, the Scordelis-Lo roof, and Mor- 
ay’s spherical shell problem. The assumed-stress hybrid 4-node quadrilateral shell element 
derived in this paper will be referred to as AQD4 (Assumed-stress Quadrilateral Direct 
4-node element) in the following discussion for convenience. The results for the present 
element are compared with the results using the QUAD4 element of the MSC/NASTRAN 
from reference [23], the Q4S element from reference [9], the AQ element from reference [10], 
the ES1/EX47, ES5/E410, and ES4/EX43 elements of the NASA Langley CSM Testbed’ 
and the AQR8 element of reference [11]. A brief description of these elements is presented 
in the appendix for completeness. The dimensions and properties for the test problems 
are chosen in consistent units. 

Patch Test 

As the first test of the accuracy of the element, the patch test problem suggested in 
reference [23] is solved. The patch is shown is Figure 3. Elements are of arbitrary shape 
patched together to form a rectangular exterior boundary. As such, boundary conditions 
corresponding to constant membrane strains and constant bending curvatures are easy to 
apply. The applied displacement boundary conditions and the theoretical solutions are also 
shown in Figure 3. The ability of the element to reproduce constant states of strains is an 
essential requirement for achieving convergence to the correct solution as the finite element 
mesh is refined. This requirement is observed by considering an individual element within 
a mesh with a complicated stress field. As the mesh is refined, the stresses within the 
elements tend towards a uniform value. Therefore, elements that cannot produce a state 
of constant strains should not be trusted to converge to the correct solution as the mesh is 
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refined [6]. The present element (AQD4) passes both the membrane and the bending patch 
tests with no error. The recovered strains and stresses are both exact. 

Straight Cantilever Beam . , 

As a second test, the straight cantilever beam problem suggested in reference [23] is solved 
for the three discretizations (six elements) shown in Figure 4. The constant and linearly 
varying strains and curvatures are produced by applying loads at the free end o e 
beam to test the ability of the element to recover these states of deformations, ihe 
theoretical results for extension, in-plane shear, out-of-plane shear, and m-plane moment 
are simply calculated from the elementary beam theory including shear deformations. 
The theoretical result for the twist is .03406, according to Timoshenko and Goodier s 
Theory of Elasticity [24]. Reference [23] quotes the answer to be .03208. Analysis with 
three successively refined meshes converged to .03385 which is much closer to the the 
theory of elasticity solution than to that of reference [23]. Therefore, the solution from 
the theory of elasticity is taken herein for normalization purposes. Normalized results for 
the present element along with the results for other elements are shown in Table 1 These 
results indicate that all elements perform well for the rectangular mesh. However, for the 
trapezoidal and parallelogram meshes which contain considerable amount of distortion, 
only the Q4S, the AQR8, and the present element (AQD4) perform well. It is noted that 
the results (except for the twist end load) are only slightly affected due to the fact that in 
this paper the stress field is expressed in the Cartesian-coordinate system. This means that 
for higher-order displacement functions it becomes less important to describe the stresses 
in the natural-coordinate system. The present element produces an error of less than 3.5% 
for all meshes and loads which indicates the insensitivity of the present element to mesh 
distortion. The present element also gives very good results for the parallelogram mesh 
with a twist end load, while the AQR8 produces an error of 15.9%. 

Curved Cantilever Beam 

Next, the curved cantilever beam problem shown in Figure 5 is solved. The beam is forme 
by a 90° circular arc. In-plane or out-of-plane loads are applied at the free end to produce 
in-plane and out-of-plane states of deformations, respectively. The theoretical solutions 
are taken to be those quoted in reference [23]. The normalized results from the present 
element are tabulated in Table 2. Results using other elements are also shown in Table 2 
for comparison. For this problem, the mesh is distorted only slightly and the results for all 
elements are good. However, the AQR8 and the present element (AQD4) perform better 
than the other elements. Once again, it is noted that the results are only slightly affected 
due to the fact that in this paper the stress field is expressed in the Cartesian-coordinate 

system. 

Cook’s Tapered and Swept Panel 

The tapered and swept panel with one edge clamped and the other edge loaded by a 
distributed shear force is analyzed next (see Fig. 6). This problem was used by Cook 
and many other researchers to test the sensitivities of finite elements due to geometric 
distortions. The panel was analyzed by a coarse 2x2 mesh and a finer 4x4 mesh. The 
reference solution for the vertical displacement at point C is taken to be 23.90 as quote 
in reference [10]. The normalized results for the present element along with the results for 
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other elements are shown in Table 3. In this problem, the mesh is distorted only slightly 
and all elements produce reasonable results. The AQR8 and the present element (AQD4) 
however, produce results that are closest to the reference solution. 

Scordelis-Lo Roof 

The Scordelis-Lo roof is shown in Figure 7. This structure is a singly-curved shell problem 
m which both the membrane and bending contributions to the deformation are significant. 
The result reported in most papers is the vertical displacement at the midpoint of the 
free-edge. The theoretical value for this displacement is quoted in reference [25] to be 
0.3086, but the normalization value quoted in reference [23] is 0.3024. The latter value is 
also used herein for normalization purposes. Because of symmetry, only one quadrant of 
the problem is modeled. The mesh on one quadrant is chosen to be NxN for N=2,4,6,8,10 
(N=number of elements along each edge) to show the convergence of the solutions for the 
present element. The results of the normalized displacement at the midside of the free-edge 
are shown in Table 4. For this problem, the mesh is made of uniform rectangular-shaped 
elements and all the elements in the table perform well. It is observed that convergence rate 
to the reference solution for the present element is roughly the same as the other elements 
and the addition of the rotational degrees of freedom does not affect the convergence rate 
of the present element for this problem. 

Morley’s Spherical Shell 

As a final test of the present element, the pinched hemispherical shell problem shown in 
Figure 8 is analyzed. The equator of the shell is chosen to be a free edge so that the 
problem represents a hemisphere loaded at four points. The load is alternating in sign at 
90 intervals, and an 18° hole is present at the top of the hemisphere to avoid needing to 
model a pole. This structure is a doubly-curved shell problem and both membrane and 
bending contributions to the deformation state are significant. Because of symmetry, only 
one quadrant of the problem is modeled. The mesh on one quadrant is chosen to be NxN 
for N=2,4,6,8,10,12 (N=number of elements along each edge) to show the convergence of 
the solutions for the present element. The results of the normalized displacements at the 
oad points are shown in Table 5. The theoretical displacement for normalization purposes 

; S A ^V° be -° 940 fr ° m reference [23]- It is seen that the AQR8 and the present element 
(ACJD4) converge to the correct solution more slowly. In fact, a 14x14 mesh produced a 
normalized result of .952 and a 16x16 mesh produced a normalized result of .972 for both 
AQR8 and AQD4. The slower convergence of these elements for this problem is attributed 
to the fact that nearly all of the strain energy is bending energy even though the membrane 
stiffness is much larger than the bending stiffness. Consequently, any small amount of 
membrane-bending coupling strongly affects the stiffness of the shell. This membrane- 
bending coupling comes about by the coupling between the normal or “drilling” rotation 
and the bending rotations by the changes in slope at element intersections [9]. The incorrect 
geometry representation causes the slow convergence for these elements. This behavior 
shows both the importance and the need for more accurate geometry representations of 
shell problems. In reference [9], MacNeal also concluded that the Q4S element converges 
slower than the QUAD4 element for this problem for the reasons discussed. However the 
Q4S does converge faster than the AQR8 and AQD4 elements for this problem. 
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Conclusions 

A simple 4-node quadrilateral shell element with 24 degrees of freedom has been developed 
which alleviates most of the deficiencies associated with such elements. The element is 
based on the assumed-stress hybrid formulation and uses the principle of minimum com- 
plementary energy. The membrane part of the element has 12 degrees of freedom and 
includes the drilling (in-plane rotational) degrees of freedom at the nodes. The bending 
part of the element also has 12 degrees of freedom. The bending part is of class C° and 
takes into account the effects of transverse shear deformations. Both in-plane and out- 
of plane displacements are assumed to have quadratic variations along the edges of the 
element, while both in-plane and out-of-plane rotations are assumed to vary linearly. A 9- 
parameter stress field is assumed for the membrane part and a 13-parameter stress field is 
assumed for the bending part. The assumed stress fields satisfy the equilibrium equations. 
The formulation of the element is simple and straightforward. The element formulation 
is derived directly for a 4-node element. This approach is in contrast to 4-node elements 
with rotational degrees of freedom which are derived from “internal” 8-node isoparamet- 
ric elements by eliminating the midside degrees of freedom in favor of rotational degrees 
of freedom at the corner nodes. This method therefore, bypasses the formation of the 
stiffness matrix for an 8-node element and the subsequent transformation of this stiffness 
matrix to that of a 4-node element, resulting in savings of computer time. Although, the 
stiffness matrix derivation is based directly on a 4-node element, most of the other element 
matrices such as the mass matrix still are derived based on an “internal” 8-node element 
which makes the derivation and implementation of the element somewhat awkward. 

The element has been demonstrated to be accurate, pass both membrane and bending 
patch tests, is nearly insensitive to mesh distortion, does not “lock”, and has no spurious 
modes. The element also has the desirable property of being invariant with respect to 
node numbering. The fact that the stresses are expressed in a Cartesian-coordinate system 
affects the results only slightly for the moderately distorted elements in the test problems 
considered. This behavior indicates that it becomes less important to describe the stresses 
in the natural coordinate system when higher-order displacement functions are present. 
Additional savings are accrued in this method by expressing the stresses in a Cartesian- 
coordinate system because, the tensor transformations of tensorial stresses in the natural- 
coordinate system to physical stresses are not performed. 

The results obtained herein are very encouraging and warrant further research to make the 
derivation of all element matrices more direct and to extend the formulation to stability 
analysis, dynamic analysis, and nonlinear analysis. 
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Appendix 

The following is a brief description of the elements (except the present element) used in 
Tables 1-5 for comparison with the present element. 

The QUAD4 MSC/NASTRAN element is a 4-node isoparametric shell element with selec- 
tive reduced-order integration. The transverse shear uses a string-net approximation and 
augmented shear flexibility [26]. This element was developed by MacNeal and is available 
in the MSC/NASTRAN finite element code. 

The Q4S element is a 4-node shell element in which the membrane part is formulated 
internally as an 8-node isoparametric element with selective reduced-order integration and 
later reduced to a 4-node element by eliminating the midside degrees of freedom in favor 
of rotational degrees of freedom at the corner nodes. This element was developed by 
MacNeal[9]. The bending part of the Q4S is the same as that of the QUAD4[9]. 

The AQ element is a 4-node assumed-stress hybrid/mixed membrane element which is 
formulated internally as an 8-node isoparametric membrane element and later reduced to 
a 4-node membrane element by eliminating the midside degrees of freedom in favor of 
rotational degrees of freedom at the corner nodes. This element was formulated by Yunus 
et al.[10]. The only result reported in reference [10] for the cantilever beam problem in 
Table 1 using the AQ element is for the mesh with trapezoidal-shaped elements with a 
unit in-plane end moment. This result is reported to be .85. The difference in the results 
between the AQ membrane element and membrane part of the the AQR8 element described 
earlier is in the selection of the assumed-stress functions. 

The ES1/EX47 element is a 4-node C° isoparametric assumed natural-coordinate strain 
(ANS) shell element developed by Park and Stanley (see, references [27-28]) and imple- 
mented in the CSM Testbed Software System[21] by Stanley using the generic element 
processor template[22]. This element is not invariant and does not pass the patch test. 
This element does not include the drilling degrees of freedom in the formulation. 

The ES5/E410 element is a 4-node C 1 shell element which was originally implemented in 
the STAGS finite element code and later in the CSM Testbed by Rankin[29]. This element 
includes the rotational degrees of freedom in the formulation and uses cubic interpolation 
for all the displacement fields. This element is not invariant and does not pass the patch 
test. 

The ES4/EX43 element is a simple 4-node ( 7 ° isoparametric assumed-stress hybrid/mixed 
shell element implemented in the CSM Testbed by this author[30]. This element passes 
the patch test and is invariant with respect to node numbering. This element does not 
include the drilling degrees of freedom in the formulation and uses linear interpolation for 
all displacement and rotation fields. 

The AQR8 element is a 4-node shell element which is formulated internally as an 8-node 
isoparametric assumed-stress hybrid/mixed element and later reduced to a 4-node element 
by eliminating the midside degrees of freedom in favor of rotational degrees of freedom 
at the corner nodes[ll]. This element was also developed and implemented in the CSM 
Testbed by this author. 
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Table 1. Normalized tip displacements in direction of 
loads for straight cantilever beam. 


Tip Loading 
Direction 

QUAD4 

MSC/ 

NASTRAN 

ESI/ 

EX47t 

ES5/ 

E410t 

ES4/ 

EX43 

AQR8t 

AQD4 

(present) 


(a) rectangular shape elements 



Extension 

.995 

.995 

.994 

.996 

.998 

.998 

In-plane Shear 

.904* 

.904 

.915 

.993 

.993 

.993 

Out-of-Plane Shear 

.986 

.980 

.986 

.981 

.981 

.981 

Twist 

.941** 

.856 

.680 

1.023 

1.011 

1.011 

End Moment 

— 

.910 

.914 

1.000 

1.000 

1.000 

(b) trapezoidal shape elements 



Extension 

.996 

.761 

.991 

.999 

.998 

.998 

In-plane Shear 

.071* 

.305 

.813 

.052 

.986 

.986 

Out-of-Plane Shear 

.968 

.763 

# 

.075 

.965 

.965 

Twist 

.951** 

.843 

# 

1.034 

1.029 

1.009 

End Moment 

— 

.505 

.822 

.102 

.996 

.995 

(c) parallelogram shape elements 



Extension 

.996 

.966 

.989 

.999 

.998 

.998 

In-plane Shear 

.080* 

.324 

.794 

.632 

.977 

.972 

Out-of-Plane Shear 

.977 

.939 

.991 

.634 

.980 

.980 

Twist 

.945** 

.798 

.677 

1.166 

1.159 

1.010 

End Moment 

— 

.315 

.806 

.781 

.989 

.986 


f These elements are not invariant and do not pass the patch test. 

J Assumed-stresses are in the natural coordinates and do not, in general, satisfy the 
equilibrium equations (see Ref. [11]). 

* The results from MacNeal’s Q4S element for in-plane shear load are reported in refer- 
ence [9] to be .993, .988, and .986 for the meshes (a), (b), and (c) in Fig. 4 respectively. 

** These results for twist were normalized with .03028 in reference £3]. Herein, all 
the other results for twist are normalized using .03046 according to Timoshenko and 
Goodier’s Theory of Elasticity [24]. 

The element produces a singular stiffness matrix for this mesh. 
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Table 2. Normalized tip displacements in direction of 
loads for curved cantilever beam. 


f— — “ — 

Tip Loading 

QUAD4 






Direction 

MSC/ 

NASTRAN 

ESI/ 

EX47 

ES5/ 

E410 

ES4/ 

EX43 

AQR8 

AQD4 

(present) 

In-plane Shear 

.833 

.929 

.938 

.888 

.997 

.996 

Out-of-Plane Shear 

.951 

.935 

.887 

.925 

.956 

.956 


Table 3. Normalized vertical deflection at point C 
for the tapered and swept panel. 


Mesh 

AQ 

ESI/ 

EX47 

ES5/ 

E410 

ES4/ 

EX43 

AQR8 

AQD4 

(present) 

2x2 

.914 

.880 

.873 

.882 

.930 

.926 

4x4 

.973 

.953 

.953 

.962 

.979 

.979 


Table 4. Normalized displacements at the midpoint of 
the free-edge for Scordelis-Lo roof. 


Mesh 

=■■ =i 

QUAD4 

MSC/ 

NASTRAN 

ESI/ 

EX47 

ES5/ 

E410 

ES4/ 

EX43 

AQR8 

AQD4 

(present) 



1.387 

1.384 

1.459 

1.218 

1.218 


1.050 

1.039 

1.049 

1.068 

1.021 

1.021 



1.011 

1.015 

1.028 

1.006 

1.006 



1.005 

1.005 

1.017 

1.003 

1.003 



1.003 

1.001 

1.011 

1.001 

1.001 
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Table 5. Normalized displacements at load points 
for hemispherical shell problem. 


Mesh 

QUAD4 

MSC/ 

NASTRAN 

ESI/ 

EX47 

ES5/ 

E410 

ES4/ 

EX43 

AQR8 

AQD4 

(present) 

2x2 

.972 

.968 

.338 

1.032 

.382 

.381 

4x4 

1.024 

1.018 

.519 

1.093 

.227 

.226 

6x6 

1.013 

1.001 

.841 

1.060 

.432 

.432 

8x8 

1.005 

.995 

.949 

1.040 

.681 

.680 

10x10 

1.001 

.993 

.978 

1.027 

.835 

.835 

12x12 

.998 

.992 

.988 

1.020 

.914 

.914 


f The drilling degrees of freedom for these elements were not suppressed in this problem. 
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y 




Location of nodes: 


node 

X 

y 

i 

.04 

.02 

2 

.18 

.03 

3 

.16 

.08 

4 

.08 

.08 


Applied 

displacements: 


(a) 

Membrane patch test 
Boundary conditions: 

u = 10- 3 (x + y/2) 
v = 10 _3 (a:/2 +y) 


Theoretical solution: 

e xx = tyy = 7 xy = 10 

<J xx — CTyy — 1333., <T x y — 400 

00 

Bending patch test 
Boundary conditions: 

w = -I0~ 3 (x 2 + xy + y 2 )/2 
0 X = — 10 _3 (a;/2 ±y) 

9 y = -10- 3 (*±y/2) 


Theoretical solution: 



Bending moments per unit length: 

M x = M y = 1.111 xio- 7 , M xy = 3.333X10- 8 

Surface stresses: 

(J xx = <Tyy = ±0.667, <r I5 , = ±0.200 

Figure 3. The patch test problem. a=0.24, b=0.12, t=0.001, E=10 6 , u=0.25. 
(Consistent units are used for various properties.) 
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s 1 1 — 1 — 1 


a) rectangular shape elements 



45 A 

X 45 

$ 

\ Z \ / 

— 1 

T~ 

b) trapezoidal shape elements 



A 45 


S — 

=Z Z A-y 

1 


c) parallelogram shape elements 


Theoretical solutions: 


Tip load direction 

Displacement in direction of load 

Extension 

1 

o 

1—1 

X 

00 

In-plane shear 

.1081 

Out-of-plane shear 

.4321 

Twist 

.03406 

In-plane moment 

.009 


Figure 4. 


Straight cantilever beam problem. Length=6., height=0.2, 
v=0.3, mesh=6xl. Loading: unit forces at the free end. 


depth=0.1 


, E=10 7 , 
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Theoretical solutions: 

Tip load direction 

Displacement in direction of load 

In-plane shear 

.08734 

Out-of-plane shear 

.5022 


Figure 5. The curved cantilever beam problem. Inner radius=4.12, outer radius=4.32, 
depth=0.1, E=10 7 , i^=0.25, mesh=6xl. Loading: unit forces at the free end. 
(Consistent units are used for various properties.) 
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(a) 



(b) 

Figure 6. The tapered and swept panel problem. Thickness=l., E=l., i/=l/3, 

mesh=NxN. Loading: unit in-plane shear force distributed on the free edge. 
Reference solution: vertical displacement at C=23.90 from reference [10]. 
(a) 2x2 mesh, (b) 4x4 mesh. 

(Consistent units are used for various properties.) 


22 




Figure 7. The Scordelis-Lo roof problem. Length=50., radius— 25., thickness=0.25, 
E=4.32xl0 8 , ^=0.3, mesh=:NxN. Loading: 90. per unit area in vertical direc- 
tion, i.e., gravity load; u x —u z — 0 on curved edges. Reference solution: vertical 
displacement at midpoint of free-edge=0.3024 from reference [23]. 

(Consistent units are used for various properties.) 
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z 



Figure 8. The spherical shell problem. Radius=10., thickness =0.04, E=6.825xl0 7 , 
i/=0. 3, mesh=NxN. Loading: concentrated forces as shown. Reference 

solution: radial displacement at the load points=0.0940 from reference [23]. 
(Consistent units are used for various properties.) 
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